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An  Aiigmented  Lagrangian  Formulation  for  the  Finite 
Element  Solution  of  Contact  Problems 

Jostph  A.  Landers^  and  Robert  L.  Thylor^ 

Department  of  Civil  Engineering, 

University  of  California.  Berkeley. 


1.  Introduction 

While  solutions  to  contact  problems  have  been  available  for  years  [l],  tech* 
niques  applicable  to  a  wide  range  of  problems  with  a  minimiim  amount  of  user 
intervention  have  been  quite  limited.  This  is  true  not  only  of  static  contact 
problems,  but  also  in  dynamic  impact  situations  as  welL 

In  this  report,  a  study  is  first  made  of  various  frictionless,  small  deforma* 
tion  contact  algorithms.  The  advantages  and  disadvantages  of  each  are  dis¬ 
cussed.  Next,  the  augmented  Lagrangian  solution  method  is  presented  in  the 
context  of  finite  element  contact  problems.  This  discussion  covers  both  the 

a 

static  and  dynanuc  algorithms.  The  generality  of  the  augmented  Lagrangian 
method  is  shown  by  examining  how  the  other  methods  may  be  derived  from  it. 
Several  examples  are  presented,  and  comparisons  are  made  to  work  done  by 
previous  researchers.  Finally,  recommendations  are  made  for  future  studies  in 
contact  problems. 
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2.  Contact  Solution  AlforiUmia 

Contact  solution  algorithms  are  generally  based  upon  either  the  classical 
Lagrange  procedure  or  the  penalty  function  method.  In  the  context  of  linear¬ 
ized  elasticity,  finite  element  solutions  can  be  viewed  as  the  minimization  of  the 
potential  energy  functional  presented  in  equation  (1). 

n(a)  =  f  W(jj)  dv  —  f  pb.jidw  -  f  t>jLda  (l) 

Where  W  is  the  stored  energy  function,  pb  is  the  body  force  vector,  t  is  the  sur¬ 
face  traction  specified  on  the  boundary  and  ji  is  the  displacement  field.  The 
integrals  are  taken  over  both  bodies  and  their  boundaries.  Application  of  stan¬ 
dard  finite  element  procedures  leads  to  equation  (2).  the  discrete  form  of  the 
functional. 

n(u)  =  )iu‘Ku-u‘R  (2) 

Where  K  is  the  system  stiffness  matrix,  f  is  the  vector  of  nodal  forces,  and  u  is 
the  displacement  vector  of  nodes  in  the  mesh.  In  the  case  of  finite  element  con¬ 
tact  problems,  a  further  requirement  states  that  penetration  of  the  two  bodies 
does  not  take  place. 

(Sl^i  U  =  0  (3) 

Where  ul^^  and  are  the  positions  in  either  deformed  body,  taken  over  the 
contact  interface,  and  n  is  the  normal  between  the  bodies. 

Application  of  equation  (2)  subject  to  the  constraint  of  equation  (3)  is 
equivalent  to  an  optimization  problem,  which  may  be  stated  as: 

minimize:  n(u)  (4.a) 

with  the  constraint:  g^u  =  0  (4.b) 

The  purpose  of  this  report  is  to  discuss  solution  schemes  for  the  expressions 

listed  in  equations  (4). 


2.1.  The  Penalty  Method 

Optimization  problems  that  involve  constraints,  such  as  those  presented  in 
equation  (4.b},  are  frequently  converted  to  unconstrained  problems  by  using 
these  constraint  equations  to  establish  relationships  between  the  unknowns  [2]. 
By  writing  equation  (2)  with  the  addition  of  a  fictitious  energy  term,  a  new  func¬ 
tional  is  defined. 

n'(u)  =  n(u)  +  ^u‘gg*u ,  tuhere;  p  >  0  (5) 

The  last  term  contains  the  constraint,  g*u,  times  a  penalty  term,  p.  Because 
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solution  to  finite  element  problems  involve  the  minimization  of  equation  (5),  the 
answer  will  tend  to  have  a  small  value  for  g*u .  Taking  the  first  variation  of  (5) 
equation  (6)  is  found. 

dtt*Ku  -  du^R -f  =  0  (6) 

By  the  fundamental  lemma  of  variational  equations,  this  becomes: 

+  (7) 

Several  comments  are  applicable  to  the  penalty  method.  First,  the  correct 
choice  of  the  penalty  parameter,  /x .  is  the  essence  of  the  algorithm.  Since  the 
penalty  approach  satisfies  the  contact  conditions  only  approximately,  valid 
results  from  equation  (4)  are  highly  dependent  on  the  value  of  y. . 


Flgxire  1:  The  effect  of  the  penalty  term,  one  dimensional  case. 


In  Figure  1,  a  plot  of  the  two  functionals  is  shown.  The  addition  of  the  penalty 
parameter  to  equation  (2),  as  stated  in  equation  (5),  can  be  seen  as  rendering 
the  functional  convex  [3].  This  is  true  even  if  n(u)  is  not  already  convex,  for 
points  near  the  solution.  Also  note  that  in  Figure  2,  the  minimum  of  the  func* 
tional  with  the  penalty  term  will  approach  the  minimum  of  the  functional 
without  the  penalty  term  only  as  /i-**. 

Issues  such  as  bring  forward  a  major  problem  with  the  penalty 

method:  the  matrices  used  in  the  solution  of  the  finite  element  problems  may 
exhibit  poor  numerical  conditioning.  While  large  penalty  values  insure  a  tighter 
tolerance  for  the  constraint,  these  values  make  an  accurate  solution  of  the 
linear  system  of  equations  much  more  difficult. 
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Figure  2:  The  effect  of  increasing  the  penalty  parameter. 


The  following  problem  will  serve  as  a  basis  of  comparison  for  the  various 
contact  solution  algorithms.  This  system  consists  of  two  springs  connected  in 
series.  One  end  is  fixed,  while  at  the  other,  a  concentrated  force  is  applied.  The 
constraint  requires  that  the  displacements  of  nodes  1  and  2  remain  equal.  Fig* 
ure  3  illustrates  the  system  and  external  loading.  By  applying  the  penalty 
method,  the  matrices  corresponding  to  equation  (7)  are: 


2k  -k 

0 

1 

K  = 

-fc  k 

R  = 

9 

g  = 

-1 

(8) 


This  leads  to  a  linear  system  of  equations; 


2k  +  fi 

“(*  +  m) 

“‘1  _  ° 

-(*  +  fi) 

k  IJL 

“zj  "  19, 

This  system  has  the  solution: 


i^: 


i 

k 

2_  2k  +  u 
k  k  +  fj, 


(9) 


(10) 


The  exact  solution  is  shown  in  equation  (11). 


SL 

“i 

k 

1^2 

i 

k. 

(11) 
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Figure  3:  A  sample  finite  element  system. 


Note  that  the  condition  number  of  the  stiffness  matrix  can  easily  be  found  to 
be: 


cond  (K)  “  —  3A;  2u  •<*  ^ SA;^  4-  -f 

3k  +Zfj,-  Vsib*  +  8fc/x  + 


As  Figures  4  and  5  illustrate,  increasing  penalty  values  lead  to  more  accu¬ 
rate  solutions,  but  they  also  lead  to  more  poorly  conditioned  matrices.  .Because 
of  the  extreme  importance  of  /i ,  a  natural  question  to  ask  is:  how  can  an 
optimal  value  for  the  penalty  parameter  fi  be  found? 

In  order  to  find  an  answer  to  this  question,  the  errors  involved  in  solution 
of  the  finite  element  contact  problem  must  be  examined  in  some  detail.  Two 
sources  of  error  effect  the  accuracy  of  the  penalty  method:  the  large  perturba¬ 
tion  that  results  from  a  small  penalty  parameter,  and  the  rounding  errors  due 
to  large  penalty  values  in  finite  precision  arithmetic  on  the  computer  [4].  A 
study  of  the  perturbation  errors  is  done  first. 

Applying  the  Sherman-Morrison  formula  to  equation  (7),  the  following 
expression  for  the  displacements  by  the  penalty  method  is  obtained. 


The  exact  solution  can  be  found  by  letting 


R 


(14) 
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For  some  constant  a,  this  expression  can  be  written  as; 


In  practice,  for  a  reasonably  large  jj. ,  a  has  a  value  near  unity. 


Log  (  Penalty  Value  ) 


Figure  5:  The  effect  of  fi  on  the  condition  number. 


A  standard  approximation  for  rounding  errors  is  given  in  [5]  as  an  expres¬ 
sion  of  the  form; 

S  b2-‘cond(K)  (21) 

where  represents  the  result  computed  using  finite  precision  arithmetic, 
identifies  the  result  computed  if  exact  arithmetic  is  used.  The  symbol  b 
represents  some  constant,  while  t  represents  the  number  of  binary  digits  of 
floating  point  precision.  On  the  DEC  VAX*  series  of  computers,  using  double 
precision  instructions,  t  has  a  value  of  56.  An  approximation  for  the  condition 
number  of  the  system  matrix  K  is; 

cond(K)  -7^  (22) 

Apjjn 

Here,  K^ri  is  smallest  diagonal  stiffness  matrix  coefficient  that  is  in  the  contact 
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region.  Hence,  the  total  error  is  represented  by  summing  both,  as  in  equation 
(23). 


Differentiating  this  expression  and  setting  it  equal  to  zero,  the  minimum  is; 


Assuming  that  the  following  approximation  holds: 

Ig'ir'gl  a  -4—  (25) 

'  *  AwriB 

equation  (24)  then  becomes: 


c  ATStof 

2-< 


(26) 


for  some  constant  c.  In  practice,  c  has  a  value  near  unity  when  /x  is  reasonably 
large.  Equation  (26)  represents  a  recommendation  for  the  penalty  value,  within 
the  limitations  presented  in  its  derivation. 

To  summarize,  penalty  methods  have  several  advantages.  Despite  the  con* 
ditioning  and  accuracy  problems,  the  resulting  solution  requires  no  extra 
efforts  beyond  adding  the  contact  conqjonents  to  the  stiffness  matrix.  Further* 
more,  these  methods  are  easy  to  implement  and  have  a  certain  intuitive  appeal 
to  engineers.  Finally,  the  contact  force  conjugate  to  the  contact  constraint, 
although  not  explicitly  found  by  this  method,  can  be  computed  using  an  extra 
procedure. 


2.2.  The  Lagrange  Multiplier  Method 

The  classical  Lagrange  methods  were  among  the  first  methods  used  to  solve 
contact  problems.  UTith  this  method,  the  functional  of  equation  (2)  is  amended 
with  an  additional  term.  This  defines  the  new  system  listed  in  equation  (27). 

n'(u,X)  =  n(u)  +  A‘g*u  (27) 

Taking  the  variation  of  (27)  yields  the  following  two  equations: 

du'Ku  -  (5u*R  +  6u*gX  =  0  (28. a) 

dX‘g*«=0  (28.b) 

By  the  fundamental  lemma  of  variational  equations,  this  leads  to  the  system; 

I*  *tl“l  _  l"l 


vnm 


! 


W 

», 

s 
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The  set  of  equations  expressed  by  equation  (29)  may  be  indefinite  but 
non-singular  [2].  Unlike  the  standard  penalty  method,  the  dimensions  of  the 
matrices  increase  with  this  approach.  Additionally,  as  illustrated  by  equation 
(29),  the  matrix  contains  a  zero  sub-matrix  along  its  diagonal.  In  a  practical 
finite  element  system,  these  zero  diagonal  terms  are  not  automatically  available 
in  a  sub-matrix,  but  are  scattered  throughout  the  system  stiffness  matrix. 
Thus,  care  must  be  exercised  when  the  solution  is  found;  special  procedures 
should  be  used  which  effectively  re-order  the  system  of  equations. 

Despite  these  disadvantages,  classical  Lagrange  multiplier  methods  allow 
the  contact  conditions  to  be  satisfied  exactly.  Also,  the  X  value  has  a  straight¬ 
forward  physical  significance:  it  is  the  value  of  the  contact  force  conjugate  to 
the  constraint  condition. 


An  illustration  of  how  the  classical  Lagrange  method  is  used  in  contact 
problems  is  demonstrated  by  again  considering  the  two  spring  problem  Apply¬ 
ing  equation  (29),  the  linear  system  of  equations  can  easily  be  found: 


2k 

-k 

ij 

' 

0 

-k 

k 

-1 

“2 

s 

? 

1 

-1 

o] 

X 

0 

This  system  has  the  solution: 


2. 

If 


(30) 


(31) 


The  results  are  the  exact  values  for  the  displacement  and  contact  force. 


2.3.  The  Perturbed  Lagrangiaa  Method 

An  interesting  extension  to  the  classical  Lagrange  method  is  offered  in 
references  [6]  and  [7].  Here  the  functional  described  by  equation  (27)  is 
modified  with  an  additional  term* 

n’(u.X)  =  n(u)  +  X‘g*u  -  '"Aere:  >  0  (32) 

The  effect  of  the  last  term  is  to  fill  in  the  zero  sub-matrix  of  equation  (29). 
The  resulting  stiffness  matrix  no  longer  contains  the  zeros  along  the  diagonal, 
the  classical  Lagrange  functional  is  recovered  by  allowing  Superficially,  it 

might  be  observed  that  this  approach  stUl  suffers  from  one  problem  of  the  clas¬ 
sical  Lagrange  method:  an  increase  in  the  number  of  degrees  of  freedom  of  the 
system  However,  this  problem  can  be  circumvented  in  a  clever  manner  when 
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the  method  is  applied  to  finite  element  systems.  To  see  how  this  may  be  done, 
consider  the  first  variation  of  equation  (32). 


dufKu  - 


+  ^Xdu‘g 


da  =  0 


(33.a) 


-X 


+  g*u 


da -0 


(33.b) 


Vfhere  the  integral  extends  over  the  area  of  contact  between  the  two  bodies. 
Since  dX  is  arbitrary  and  fj.  is  a  given  constant,  equation  (33. b)  can  be  written 
as: 


j-f  Xio  =/g*uda 

^  <e  *• 


(34) 


In  the  framework  of  finite  element  analysis,  the  displacement  field  is  generally 
approximated  over  each  element  region  by  some  polynomial.  If  this  approxima¬ 
tion  is  linear,  the  strain  field  is  constant  function.  Thus  assuming  that  the  con¬ 
tact  force,  X,  is  constant  over  the  region  of  interest,  the  above  expression  can 
be  transformed  to: 


X=  -f-fg^uda 


(35) 


Because  of  the  linear  approximation,  the  integral  can  be  evaluated  by  means  of 
the  trapezoidal  rule.  Hence,  X,  can  be  approximated  to: 

(38) 


The  point  of  this  exercise  is  that  equation  (33. a),  after  applying  the  funda¬ 
mental  lemma  of  variational  equations,  becomes  the  expression  presented  in 
equation  (37). 

Ku  >  R  gX  =  0  (37) 

Applying  the  results  of  (36).  a  new  expression  is  found. 

Ku-R  +  Zigg'u^  0  (38) 

Thus,  the  two  equations  in  (33)  can  be  reduced  to  a  single  system  that  does  not 
increase  the  number  of  degrees  of  freedom.  Equation  (38)  is  essentially  that  of 
equation  (7),  except  for  the  following  points.  First,  an  explicit  expression  for 
the  contact  force  is  directly  available  through  the  use  of  equation  (38).  Second, 
while  a  solution  is  sought,  using  the  Lagrange  approach,  to  the  following: 

minimize:  n(u)  (39.a) 

with  the  constraint:  g*u  =  0  (39.b) 

What  instead  is  obtained  by  using  the  perturbation  is: 


minimize :  n(a)  (40.a) 

with  the  constraint:  =  0  (40.b) 

Clearly,  one  can  see  that  the  solutions  are  equivalent  only  as  but  as  was 
seen  in  the  standard  penalty  approach,  increasing  fj.  leads  to  numerical 
difficulties. 


Again,  the  two  spring  problem  is  presented  within  the  context  of  the  per¬ 
turbed  Lagrangian  method  in  order  to  illustrate  how  it  may  be  applied  to  finite 
element  systems.  The  system  of  linear  equations  is: 


2k  +fj. 

-(X  +  m) 

“il 

fo 

”(*:  +  m) 

fc  +  /i 

Ug] 

~  P, 

(41) 


This  has  the  solution: 


1 

k 

u 

“2 

i 

2k  +  u 

A - q 

[k  fi 

k 

k  +  n 

(42) 


Note  that  as  the  exact  solution  is  found. 


2.4.  The  Augmented  Lagrangian  Method 

Until  now,  the  methods  considered  have  had  some  considerable  implemen¬ 
tation  drawbacks.  The  penalty  and  perturbed  Lagrangian  methods  had  solu¬ 
tions  which  are  highly  dependent  on  the  penalty  parameter  fj, .  Although  it  is 
possible  through  the  use  of  an  expression  like  that  presented  in  equation  (26) 
to  get  a  reasonably  good  value  for  fi  based  on  conditioning  considerations,  such 
a  value  says  nothing  about  how  accurate  a  given  solution  might  be.  Moreover, 
the  classical  Lagrange  method,  while  it  promises  the  exact  results,  requires  that 
the  resulting  larger  system  of  equations  effectively  be  re-ordered  during  the 
solution  of  the  system.  One  might  ask  if  it  is  possible  to  develop  a  method 
which,  while  retaining  a  relatively  easy  implementation,  minimizes  the  disadvan¬ 
tages  of  the  penalty  and  Lagrange  methods.  To  this  end.  consider  the  following 
functional  in  equation  (43). 

n'(u.X)  =  n(u)  +  X‘g*u  +  ^u‘gg*u.  where  ifj.  >  0  (43) 

Equation  (43)  can  be  viewed  as  a  classical  Lagrangian  method: 

minimize:  n(u)  +  ^u*gg*u 
with  the  constraint:  g*u=  0 


(44.a) 

(44.b) 


because  the  addition  of  the  penalty  term  to  equation  (44. a)  does  not  change  the 
optimal  values  of  n  or  the  Lagrange  multipliers  X .  Also,  note  that  the  system 
can  also  be  seen  from  the  penalty  viewpoint: 

minimuie :  n(u} -f  X'f'u  (45. a) 

with  the  constraint:  g^u  =  0  (45. b) 

Observe  that,  if  the  correct  values  of  the  Lagrange  multipliers  X*  ere  used, 
equations  (45)  can  be  written: 

mintmiee:  n(u)  +  (X')*g*tt  (48.a) 

with  the  constraint:  ^u-0  (48. b) 

Or,  by  applying  the  penalty  approach: 

n*(u.x)  =  n(u)  +  (x*)*g‘u  +  ^u‘gg*u  (47) 

Taking  the  first  variation  of  equation  (47)  and  noting  that  since  X*  is  the  exact 
value,  the  result  must  be  equal  to  zero: 

du*Ku  “  du^R  +  du'gX*  + /idu'gg^u  =  0  (48. a) 

(dX*)^g*u=0  (48. b) 

Since  du*  and  (<5X*)*  are  arbitrary,  equations  (48)  have  the  solution: 

Ku  -  R  +  gX*  =  0  (49.a) 

g*u=0  (49.b) 

The  formulas  (49)  represents  a  set  of  equations  for  an  exact  penalty  method. 
For  details  on  exact  penalty  methods,  see  [3]. 

Since  the  correct  values  of  X  are  not  known  in  advance,  a  procedure  must 
be  developed  to  find  them,  if  the  functional  of  equation  (43)  is  to  be  of  any  use. 
Until  now,  the  algorithms  developed  did  not  require  iteration  if  the  contact  area 
between  the  two  bodies  remained  constant.  Tith  the  augmented  Lagrangian 
method,  an  iteration  is  required  to  conq>ute  the  correct  values  of  u  and  X  .  At 
first  this  procedure  might  seem  to  be  at  a  disadvantage,  since  extra  computa¬ 
tion  is  required.  However,  as  will  be  shown,  there  are  several  advantages  to  this 
approach.  The  system  can  yield  a  value  for  X  within  a  given  tolerance,  and 
hence,  u  can  be  computed  to  a  specified  value.  This  did  not  occur  with  the 
penalty  or  perturbed  Lagrangian  methods.  Second,  the  number  of  equations  of 
the  S3rstem  does  not  increase.  This  was  not  the  case  in  the  classical  Lagrangian 
method.  Third,  convergence  to  a  specified  value  does  not  require  that 
This  can  help  mitigate  the  severe  conditioning  problems  associated  with  using  a 
large  penalty  value.  Fourth,  and  most  important,  in  finite  clement  problems, 
the  contact  area  is  generally  not  known  in  advance.  This  area  must  be  found  as 
the  solution  progresses.  Hence,  an  iterative  procedure  is  always  necessary  for 


w. 

(45.a) 

(45.b) 

V 

are  used. 

(48.a) 

.S. 

^S* 

(48.  b) 
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the  finite  element  solution  by  the  penalty,  Lagrangian,  and  perturbed  Lagran- 
gian  methods  previously  presented. 

A  first  order  update  method  for  X  can  be  written  as  in  equation  (50). 

A*>i  =  X*+Mg*u  (50) 

The  motivation  for  such  an  updating  scheme  can  be  found  if  equation  (45)  is 


written  as: 

mimirnize :  n(u)  -f  X^g^u 

(51.a) 

with  the  constraint:  g*u  =  0 

(51.b) 

Viewing  this  from  a  Lagrange  standpoint: 

n(u)  +  xig‘tt+  (X*  -x*)‘  g‘u  =  0 

(52) 

A  typical  step  in  the  augmented  Lagrangian  method  involves  the  minimuation 
of: 


n(u)  +  Xig;*a  +  =  0 

So,  comparing  equations  (52)  and  (53).  a  good  approximation  is: 

-X* 


(53) 

(54) 


Rearranging  terms,  Xjb+i  is  found. 

=  X*  +  Mg*u 


(55) 


This  is  exactly  the  updating  method  proposed  by  equation  (50). 

Considering  the  two  spring  problem  described  in  Figure  3  for  the  aug¬ 
mented  Lagrangian  method,  the  following  system  of  linear  equations  results: 


(56) 


Zk  +  fi  —{k  +  fj,) 

Ul] 

C 

-(k  +  ft)  k+fi 

UgJ 

"  h 

rj"^ 

|-i 

Note  that  in  this  example.  X  is  one  dimensional  and  that,  initially  it  is  set  to 
zero. 


Xo  =  0  (57) 

So,  the  equation  to  update  X  is  simply: 

X**i  =  X*  +m(wi  -  Wz)  (58) 

Solving  the  linear  system,  and  updating  the  value  of  X  leads  to  the  following 
series  of  equations. 


^1 


-Jil. 

k  +  fj. 


Xj 


k 

k  •¥  n 


(59.a) 

(59.b) 
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> 


Xg  =  T-^1  +  r-==--  + 
^  *  +  At]  ib  +  /i 


k 

k  +  pt 


(59.c) 


Or,  in  general,  equations  (59)  can  be  written  as: 


X_  =  ^  k  ^  k  k  k 

^  k  jj.  k  +  fj,  k  +  fi  k  fi  .  k  +  fi 


This  has  the  form: 

X„  =  50(1  +  a  +  a®  +  o’  + . +  a’‘-‘) 

for  the  following  substitutions: 

^  k  +  fi  k  +  jj. 

An  infinite  geometric  series,  this  converges  to: 


X  •  - ^  1  _  -M 

1 

^  1  —  a  k  +  fi 

1-  * 

k  +  fi 

A  requirement  for  convergence  is  that: 

LI.,  ..  L 

k  L, 

1  A  W#  ,  ^  ’S.  A 

11  k  ■¥  fi 

Therefore,  it  is  seen  that  this  system  will  converge  for  ony  positive  value  of  fj, . 
At  convergence,  the  displacements  can  easily  be  found  as: 

L  1  ^  i  ^ 

-  *  _  * 

“zj  "  0  >  X’hk  ^  ~  3. 

k  k  +  fj,  k  * 

These  values  represent  the  exact  result. 

For  a  detailed  convergence  proof  of  the  augmented  Lagrangian  method  in 
the  general  case,  see  [3].  The  example  previously  presented  displays  the  power 
of  the  method  for  the  following  reasons.  First,  convergence  to  a  specified  toler* 
ance  is  found  without  requiring  Second,  the  correct  displacements  and 

contact  forces  are  available  without  the  need  for  extra  calculations.  Also  note 
that,  in  practice,  the  penalty  parameter  does  not  have  to  be  held  fixed  during 
the  iteration:  it  can  slowly  increase  as  the  solution  nears  convergence.  Third, 
the  burden  of  accuracy  does  not  depend  directly  on  /j.  but  rather  on  the 
number  of  iterations  used  to  compute  A .  Furthermore,  the  algorithm  does  not 
have  to  be  run  until  X*  is  found,  but  rather  it  needs  only  to  be  computed  until  X 
is  within  a  specified  tolerance  of  X*.  Finally,  each  iteration  leads  to  an 


improved  estimate  of  the  solution. 

In  order  to  study  further  the  convergence  properties  of  this  method,  the 
duality  viewpoint  of  equation  (43)  will  be  examined.  Recall  that  duality  means 
that  the  functional  which  results  from  the  problem; 

minimize :  n(u)  +  ^u'gg'u  (86.a) 

with  the  constraint:  g*u  =  0  (66.b) 

will  lead  to  the  system  of  equations; 

Ku  -  R  +  {X*)*g  +  Atgg*u  =  0  (67) 


when  the  Lagrange  approach  is  used.  By  invoking  a  convexity  assumption  on 
the  Hessian  of  equation  (67),  an  equivalent  dual  problem  can  be  written: 


jj(X)  =  minimum 


n(u)  +  XV^  + 


(88) 


Specifically,  the  dual  function  is  used  since  it  is  a  function  of  only  one  variable, 
X  ,  and  hence  there  is  a  corresponding  simplification  in  the  problem  statement. 
The  dual  problem  requires  that  the  Hessian  of  the  solution  of  the  solution  to 
equations  (66)  be  positive  definite;  this  is  a  convexity  assumption.  Note  that 
from  the  results  of  the  penalty  method,  and  in  fact  from  the  results  of  Figure  1, 
the  problem  is  locally  convex  near  the  solution  point  for  a  reasonably  large 
value  of  fx  .  Also  worth  noting  is  that  equation  (68)  is  not  used  in  the  solution  to 
the  finite  element  problem  since  the  minimization  requires  a  solution  of  an 
unconstrained  problem  in  terms  of  the  unknown  displacements.  Updating 
schemes  for  u  corresponding  to  equation  (50)  are  much  more  difficult  to 
develop  and  use. 

Using  f(x(X))  to  represent  the  standard  finite  element  system,  h(x(X))  to 
represent  the  constraint,  and  x(X)  as  the  vector  which  minimizes  equations 
(66),  the  dual  of  equation  (66)  can  be  written  as; 


d(X)  =  /(x(X))  +  X‘h(x(X))  +  |-h‘(x(X))h(x(X)) 


(89) 


Equation  (69)  has  the  gradient: 


Vd(X)  =  7xx(X) 


V,f(x(X))  +  V.h(x(X))X  +  MV.h(x(X))h(x(X)) 


+  h(x(X)) 


(70) 


=  h(x(X)) 

The  term  in  braces  vanishes  at  the  stationary  point,  because  it  corresponds  the 
first  variation  of  the  functional  given  by  equation  (45).  The  Hessian  of  the  dual 
is*  now  calculated  in  order  to  examine  the  convergence  properties  of  the 
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augmented  Lagrangian  method. 


V2d=Vxx{X)V,h(x(X)) 

If  the  term  in  braces  of  equation  (70)  is  described  as: 

V,L(*(X).X)  =  V.f(x(X))  +  X‘V,h(x(X))  +  MV.h(x(X))h(x(X))  =  0 
then,  upon  differentiating  this  expression  one  finds  (73). 

Vxx(X)V|,L(x(X).X)  +  V^L(x(X).X)  =  0 
Equation  (73)  can  be  further  simplified  to: 

Vxx(X)V|,L(x(X).X)  +  7,h‘(x(X))  =  0 
Finally,  Vxx(X)  can  be  found  by  solving  equation  (74). 

Vxx(X)  =  -V.h‘(x(X))V|,L-‘(x(X),X) 

Substituting  this  result  into  equation  (71).  the  Hessian  is  listed  in  (76). 
V^d(X)  =  -V,h‘(x(X))V^L-‘(x(X).X) 


=  -V,h‘(x(X)) 


yi,f(x(X))  +  AiV.h(x(X))y,h‘(x(X)) 


-1 


(7lj 

(72) 

(73) 

(74) 

(75) 

(76) 


V.h(x(X)) 


The  steepest  ascent  formula  for  maximizing  <i(X)  of  equation  (69)  can  be 
written  as: 

X**,  =X* +  aVd(X*)  (77) 

Or,  by  substituting  the  results  of  (70),  equation  (77)  becomes: 


+  ,  =  Xfc  +  oh(x(Xt)) 


(76) 


And.  from  reference  [3],  a  relationship  between  the  values  of  X  and  the  conver 
gence  ratio  is: 


-  A* 
X*  -X* 


«  r(/i)  =  max  ( 


1 


1 


—  fJ,W 


l> 


(79) 


Where  W  is  the  maximum  eigenvalue  of  V^d{k) ,  w  is  the  minimum  eigenvalue  of 
V^d(X)  ,  the  function  r(/x)  is  the  convergence  ratio  and  is  the  optimal  value  of 
/X  .  This  optimal  value  can  be  computed  to  be: 


# 


2 

W  +  w 


(80) 


And,  the  convergence  ratio  at  that  value  is: 


W-w 
W  +  w 


(81) 


Similarly,  the  Hessian  of  the  primal  function  can  be  computed  to  be; 


-Ml 


Substituting  and  simplifying,  the  Hessians  of  both  the  dual  and  primal  function¬ 
als  become; 


(X)  =  -g* [k  +  /igg*]  g 


V®p(u)  = 


g‘(K+Mgg*)  g 


-Ml 


(83) 

(84) 


Next,  consider  the  following  rank  one  update  formula: 

Mg‘(K  +  Mgg‘’“’g  =  I-(I  +  Mg*K'*g)"*  (85) 

Multiplying  both  sides  by  (I  +  /x^IC~'g)  on  the  right  leads  to  equation  (86). 

Mg*(K  +  Mfifi‘)~‘e(I  +  =  Mg*K~*g  (86) 

Suppose  e  is  an  eigenvector  of  g*K“‘g  with  an  eigenvalue  a.  Post-multiplying 
equation  (86)  by  e  and  simplifying,  it  can  be  written  as  equation  (87). 

Mg‘(K  + Mgg‘)"‘g(l  +  M2)e  =  M“e  (87) 

This  equation  has  the  form: 

fiA(l  + /u.a)e  s  fxae  (88. a) 


or; 

.  A<be  =  ae  (88. b) 

So,  e  is  also  an  eigenvector  of  A.  Thus  g*(K  +  Mgg*)**g  and  g*K'*g  have  the  same 
eigenvectors.  Using  this  information  in  equation  (87)  and  the  symbol  b  to 
represent  and  eigenvalue  of  g‘(k  +  Mgg*)~*g'  ®  relationship  between  eigenvalues 
a  and  b  can  be  written: 

6(l  +  mo)  =  a  (89.a) 

or: 

b  = 

l+MO 

By  substituting  these  values,  the  convergence  ratio  is: 


(89.b) 


Clearly,  as  the  convergence  ratio  goes  to  unity,  implying  arbitrarily 

fast  convergence.  By  application  of  the  analysis  of  the  penalty  method,  it  is 
known  that  from  equations  (13)  to  (26).  rounding  errors  also  play  a  part.  Thus, 
a  tradeoff  situation  has  developed.  Although  increasing  the  penalty  parameter 
increases  the  rate  of  convergence,  it  also  leads  to  more  poorly  conditioned  set 
of  equations.  Since  the  augmented  Lagrangian  will  converge,  within  certain  res¬ 
trictions,  for  any  value  of  the  penalty  parameter  it  is  clear  that  "best"  value  to 
use  is  the  one  which  wUl  minimize  the  rounding  errors.  Thus,  this  optimal  value 
for  the  penalty  parameter  has  exactly  the  same  form  as  equation  (26). 

As  a  further  justification  for  the  first  order  update  given  by  equation  (50), 
consider  the  application  of  Newton’s  method,  and  an  approximation  to  the  Hes¬ 
sian  of  the  dual  functional.  The  value  of  can  be  thought  of  as  minimizing 
the  second  order  expansion  of  the  dual  around  X^: 

d*(A)  =  d(X*)  +  7d‘(Xfc)(X  -X*)  +  3$(X  -X*)‘V2d(X*)(X-X*)  (91) 

Taking  the  derivative  and  setting  it  equal  to  zero: 

Vdfc(X)  =  Vd(Xfc)  +  V2d(Xfc)(X  -  X*)  =  0  (92) 

This  leads  to  the  update  formula: 

X**,  =  X*  -  [v2<i(X*)pVd(Xj  (93) 

By  equation  (70),  Vd(Xfc)  is  known,  and  from  equation  (76)  72ci(X*)  is  also 
known.  Equation  (83)  represents  the  system  of  equations  for  the  contact  prob¬ 
lem:  this  expression  could  be  substituted  into  equation  (93).  As  a  further 
simplification,  note  that  for  a  reasonably  large  /u ,  equation  (83)  is  approxi¬ 
mately  .  Hence,  equation  (93)  can  be  expressed  as: 

=  X*  +  Atg'u  (94) 

This  is  exactly  the  method  proposed  in  equation  (50). 

The  augmented  Lagrangian  procedure  can  be  interpreted  geometrically  in 
Figure  8.  Consider  the  primal  functional  given  in  equation  (43).  By  addition  of 
the  penalty  term,  the  system  is  convex  near  the  solution,  u*  ,  X*  .  For  a  given 
set  of  starting  values  for  Uq  and  Xo  the  solution  proceeds  as  follows.  The  value 
of  lit  is  found  and  is  used  to  update  the  Lagrange  parameter  X^^j  .  This  value 
corresponds  to  the  negative  slope  of  the  functional  at  the  point  ■  The  new 
value  of  Xt4.|  is  then  used  to  compute  a  new  Ut  +  i  .  This  process  may  be  contin¬ 
ued  until  the  values  Ut+n  and  Xt+„  are  within  a  specified  tolerance  to  u*  and  X*  . 


Figure  6:  Grar  hical  interpretation  of  the  augmented  Lagrangian  method. 


3.  Finite  Element  Implementation 

The  finite  element  method  may  be  used  to  construct  solution  procedures 
for  the  four  contact  schemes  discussed  in  Section  2.  In  this  section,  only  the 
frictionless,  two-dimensional,  small  deformation  kinematics  case  will  be 
presented.  Since  the  algorithms  for  the  penalty  and  classical  Lagrange 
methods  are  so  well  known,  the  discussion  of  these  methods  will  be  confined  to 
their  relationship  with  the  augmented  Lagrangian  scheme.  Details  on  the  per¬ 
turbed  Lagrangian  algorithm  are  given  in  references  [6]  and  [7]. 

The  augmented  Lagrangian  method  has  been  implemented  in  the  finite  ele¬ 
ment  program  FEAP  [9].  Details  of  the  contact  kinematics  can  be  found  in  [6] 
and  [7].  In  these  references,  the  contact  interface  is  first  split  into  a  number  of 
segments  along  an  intermediate  surface,  as  shown  in  Figure  7.  Assuming  four 
node  isoparametric  elements,  these  segments  are  unambiguously  defined  by  the 
geometry  of  the  adjacent  elements  of  the  bodies  in  contact.  Figure  6  shows  a 
typical  segment. 


Figure  8;  A  typical  contact  segment. 


lixj  -  xjll 

From  equations  (96),  tbe  displacement  can  be  described  a^s: 

u*  =  (1  -  a*)ui  + 

And,  hence,  the  gaps  in  the  contact  segments  can  be  represented  by: 
=  (u|  -  u  *).n* 

=  (W  -  u®).ii* 

The  normal  n*  is  defined  by: 
n‘  =  es  ®  t‘ 


(96.c) 

(97) 

(9B.a) 

(98.b) 

(99) 


By  applying  a  linear  interpolation  to  the  displacement  field  described  by 
equation  (97),  a  consistent  approximation  is  obtained  with  the  four  node  iso¬ 
parametric  element.  In  the  context  of  this  approximation,  the  contact  pressure 
within  each  segment  is  a  constant.  Note  that  this  approach  is  not  based  on 
nodal  penetration  of  the  adjacent  body,  but  rather  an  average  penetration 
across  a  contact  segment.  Thus,  in  the  solution  to  a  finite  element  contact 
problem,  the  constraint  is  satisfied  in  an  "average"  sense  over  a  segment 
instead  of  "point-wise"  for  a  contact  node.  By  using  the  expressions  in  equation 
(98)  for  the  interpolation  functions  g ,  in  any  of  the  equations  (5),  (27),  (32)  or 
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(43)  for  g*u,  this  averaging  can  be  incorporated  into  the  solution  methods  pre¬ 
viously  discussed. 

Consider  equation  (5),  which  represents  the  penalty  method.  Assuming  two 
bodies  in  contact,  the  functional  can  be  written  as  in  equation  (100). 


n-(u‘.u2)  =  n(u‘)  +  n(ii2)  +  ^  f-g'g  da 


(100) 


Where  the  first  two  terms  of  equation  (lOO)  are  of  the  same  form  as  equation 

(1): 

U(u.i)  =  f  Will)  dv  i  -  f  pb  •  tx  dv  I  -  f  t.udoi  (101. a) 

n( Ug)  =  rw(u)  dt(2  -  Tpb-udvj  -  f  t-udcg  (101. b) 

By  substituting  equation  (98)  into  equation  (100),  the  following  expression  is 
obtained: 


n'(u)  =  n(u)  +  ^ ^  ||u*-u*|*nj  |u®-u‘j.n 


The  kinematically  admissible  variations  are: 


V*  =  |i7,  :  ▼  ‘(Qj)-*R®  .  Wi  I  dS,  =  0 
V®  =  ^2:'r2(Q2)-R»,i,2  I  dP,  *  o| 


(102) 


(103.a) 


(103.b) 


Applying  the  Gateaux  derivative  to  equation  (102),  and  invoking  stationarity: 

D,  1  n.»7i  =  ■T^n(v‘  +  ai7,)|«,o*  +  o»7i)la.o  +  (I04.a) 

*•1  ov  a 

w 

f  {STT  ■f' 

||v  ®  +  aijgj  -  |v  *  +  ai7,j | .  n  do  =  0 

'  '  ■  0 

D,  ill. 172  =  -:^n(v®  +  ai72)|„.o‘  +  a’7z)la.o  +  (l04.b) 


IL 

2 


+  aii8 


+  ai7i 


0 


Or.  simplifying  and  combining  these  two  equations  the  following  expression  is 
obtained. 


•^^,1  n.»7i  +  D^t  n.»j2  +  /X 


da  s  0 


(105) 


Applying  the  linear  finite  element  discretization,  equation  (105)  becomes: 


+  E  p  ^  *  Ct 


t  I  i>l 


=  0 


(108) 


Where  the  first  sum  represents  the  standard  contributions  of  stillness  and  force 
to  the  finite  element  system,  and  the  second  sum  represents  the  contact  con¬ 
straint  contribution.  The  vector  c*  can  be  expressed  with  reference  to  Figure  8 
as  follows. 


Cl  a  -Ji(l  -a>)n‘ 

(107.a) 

Ca  a  +H(tt®  “  a‘n‘) 

(107.b) 

ca  a  +J$(n*  -  a'n*) 

(107.C) 

C4  =  -)6(1  -  a*)  n* 

(107.d) 

In  equation  (108),  the  sum  is  taken  over  all  of  the  elements  representing  the 
discretization  of  both  bodies  for  the  first  term,  and  the  sum  over  all  segments  in 
contact  for  the  second  term.  A  typical  contribution  of  a  contact  element  to  the 
system  stillness  matrix  can  be  described  as:  . 


Cl 

Cs 

Cs 

C4 


jCj  Cg  Cj  C4J 


(108) 


Which,  for  each  segment  in  contact,  represents  a  rank  one  update  to  the  system 
stillness  matrix. 


Some  difiUculties  are  present  with  this  kinematic  formulation,  particularly 
when  higher  order  interpolations  are  used  or  an  extension  to  three  dimensions 
is  made.  The  small  deformation,  linear  kinematics  implementation  has  a 


reasonable  geometrical  deflnition.  as  described  in  Figures  7  and  B,  for  most 
cases.  However,  the  extension  to  multiple  contact  segments  along  a  single  slide 
segment,  or  the  general  application  of  symmetric  boundary  conditions  is 
presently  unclear.  Figure  9  shows  some  problems  that  may  arise  even  for  the 
four  node  case.  Here,  the  interpolation  method  breaks  down  because  projec¬ 
tions  are  not  simple  enough  to  be  described  by  the  application  of  of  the  method 
shown  in  Figure  B. 


Figure  9:  Possible  problems  with  the  slide  line  logic. 


Furthermore,  ambiguity  arises  in  the  case  of  higher  order  interpolations  of 
the  intermediate  surface.  'While  this  problem  could  be  overcome  by  a  complex 
calculation,  it  is  not  clear  that  the  test  for  penetration  could  be  done  efficiently. 
For  example,  determining  if  a  node  is  penetrating  a  body  requires  the  use  of  a 
concave  clipping  algorithm.  For  boundaries  represented  by  polynomials  greater 
than  first  order,  the  calculation  can  be  very  expensive  [10].  Moreover,  in  order 
to  extend  this  approach  to  three  dimensions,  techniques  from  constructive 
solid  geometry  (CSG)  would  have  to  be  employed.  This  means  that  the  penetra¬ 
tion  calculations  might  generally  take  longer  than  the  solution  of  the  system  of 
equations  for  the  entire  finite  element  system  [llj- 

The  finite  element  implementation  of  the  augmented  Lagrangian  method 
can  now  be  discussed  in  more  detail.  By  applying  the  procedure  given  in  the 
last  section  to  the  functional  expressed  in  equation  (43),  a  finite  element  sys¬ 
tem  of  the  form  of  equations  (106)  can  be  derived.  This  system  is  shown  in 
equation  (I09.a)  in  matrix  notation. 


ilPu*  -  r* 


£ 
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c*{c*)‘u  +  CgXJl 


=  0 


(109.a) 
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The  equation  for  updating  X  can  be  written  as: 

X2*i  =  X| +m(c*)‘u  (lOS.b) 

Several  features  are  worth  repeating  about  equations  (l09).  First,  the 
summation  extends  over  all  of  the  elements  representing  the  discretization  of 
the  bodies  in  contact.  Second,  the  contact  stillness  contribution  is  a  rank  one 
update  of  the  global  stillness  matrix  over  each  segment  in  contact.  Third,  the 
method  is  inherently  iterative  because  the  Lagrange  parameter  is  not  known 
before  the  solution  is  attempted.  Fourth,  the  interpolation  across  a  contact 
segment  is  done  in  an  "average"  sense,  rather  than  "point-wise."  Finally,  note 
that  even  though  the  solution  is  iterative  due  to  the  nature  of  equation  (109),  it 
is  also  not  known  in  advance  which  segments  are  in  contact:  each  loading 
configuration  requires  control  logic  based  on  the  current  values  of  the  displace¬ 
ments.  Newton's  method  is  used  to  solve  equations  (109).  The  solution  pro¬ 
cedure  is  outlined  as  Algorithm  1. 

In  practice,  each  contact  segment  has  information  containing  the  current 
interpolation  functions  c,  current  Lagrange  parameter  X,  and  penalty  value  fi . 
Generally,  the  penalty  value  varies  slowly  from  its  original  value  to  something 
about  ten  to  a  hundred  times  larger  as  the  solution  converges.  This  has  the 
effect  of  tightening  the  constraint  as  the  answer  becomes  better  known  while 
keeping  the  solution  reasonably  well  conditioned  at  the  early  stages. 

The  fact  that  the  augmented  Lagrangian  method  is  a  combination  of  both 
the  penalty  and  Lagrangian  methods  has  already  been  shown.  Depending  on 
ones  viewpoint,  it  is  relatively  straightforward  to  extend  or  simplify  Algorithm  1 
to  incorporate  these  cases.  For  the  penalty  method,  the  Lagrange  parameters 
can  obviously  be  ignored,  leading  to  a  system  of  equations  analogous  to  equa¬ 
tion  (109): 

£  (K»u-R*)  +  2mc*(c*)‘u=  0  (110) 

«  « 

A  similar  procedure  holds  for  the  perturbed  Lagrangian  method,  since  equation 
(36)  expresses  a  value  for  the  Lagrange  parameter  across  the  contact  segment. 
Here,  the  step  in  Algorithm  1  which  updates  the  Lagrange  parameter  may  be 
skipped  since  this  item  is  considered  constant  across  a  given  load 
configuration. 

Step  (II.D)  in  Algorithm  1  states  that  the  exact  value  of  X  does  not  have  to 
be  found.  The  user  can  select  the  desired  precision.  Furthermore,  this  preci¬ 
sion  is.  unlike  that  of  the  penalty,  classical  Lagrangian  or  perturbed  Lagrangian 
methods,  independent  of  the  penalty  value.  This  penalty  value,  by  the  results  of 
section  2.4,  governs  the  rate  of  convergence  of  the  augmented  Lagrangian 
method.  It  is  generally  expected  that  results  to  machine  precision  can  be 
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static  Augmented  Lagrangian  Method 
I.  Set  uo  =  0.  Xq  s  0 

n.  For  k  =  1,  2,  3, ...  until  convergence  do: 

A.  BuUd  the  system  and  I^t,m 

B.  For  s  =  1  to  Stttui 

1.  Determine  e*  ,g*  =  (c*)*Ufc 

2.  Check  to  see  if  the  segment  is  in  contact: 

If  L*  >  — 1 


then  the  segment  is  not  in  contact 
else 

a.  Calculate  the  stiffness  contribution: 

b.  Calculate  the  residual  contribution:  Contact  ~  ®* 

c.  Add  the  contact  K^gtitost  •  ^ojuaet  system  l^ygttrn  • 
to  get  "Ktatal  •  Ktolot 

C.  Solve  the  system  for  the  displacement  increment: 

JQoJot  ”  R«o<ai  JQotoi 

D.  Check  for  convergence: 

if  llKtoiaj  Ufc  -  R«o«aJ  I!  <  tolerance 
then  announce  convergence  and  stop. 

E.  Update  the  displacement  field:  s  +  Au^^-i 

T.  For  s  =  1  to  Stotai  • 

Update  the  Lagrange  parameters:  X*  +  i  =  X*  + 

G.  Next  k,  go  back  to  step  II. 


Algorithm  1:  Description  of  the  static  augmented  Lagrangian  method. 


obtained  after  only  a  few  iterations. 

Since  the  algorithm  depends  on  determining  the  state  of  a  contact  segment 
in  step  (I1.B.2),  it  should  be  noted  that  the  convergence  properties  discussed  in 
section  2.4  might  not  always  apply.  This  is  because,  as  segments  change  state 
from  "not  in  contact"  to  "in  contact"  or  vice-versa,  the  tangent  matrix  becomes 
discontinuous.  For  problems  involving  practical  engineering  applications,  these 
discontinuities  last  only  during  the  initial  iterations  and  the  solution  quickly 
converges. 

The  augmented  Lagrangian  method  can  easily  be  extended  to  handle  prob¬ 
lems  involving  dynamic  contact.  Algorithm  2  presents  a  solution  of  the  problem 
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Dynamic  Augmented  Lagrangian  Hethod 
L  Initialize  °u.  °u.  °ii 
n.  For  t  s  At.  2Ai,  3Ai,  ...  until  T^tnot 

A.  Set  Uq,  Xq  =  0 

B.  For  k  s  1.  2,  3, ...  until  convergence  do: 

1.  BuUd  the  system  and 

2.  For  s  =  1  to  do: 

a.  Determine  c*  .  g*  =  (c*)‘uyE 


b.  Check  to  see  if  the  segment  is  in  contact: 


then  the  segment  is  not  in  contact  else 
L  Calculate  the  stiffness  contribution:  K^gniaet  =  Mc*(c*)*u^ 
iL  Calculate  the  residual  contribution:  Beeniwi  -  ~(Xj^'c* 
ilL  Add  the  contact  KSantact  .  Bfoiuoci  to  system  *K,y,um  • 
to  get  Ktotal  •  BloioX 

3.  Build  ^Biynamie  “  ^BtysUm  ^ 

4.  Calculate  the  residual: 

*Rd»n«mic  *‘R»olai  - ‘Ki„.t‘u- *11^,1^  0,[u*  -  *uj  -  02*6  -  Oa^U 

5.  Solve  the  system  for  the  displacement  increment:  'I^ynomic  s 

6.  Update  the  current  configuration:  u*^!  -  u^  -f  fiu 

7.  Check  for  convergence: 

if  -*K«o«az«*+i  <  folsronce]  then 

a.  announce  convergence  in  k:  *^*u  - 

b.  Update  the  velocity  and  acceleration  fields: 

‘**u  =  ao[**‘u  -  ‘uj  -  02*u  -  a3*u 

=  *U  +  Oj*U  +  07‘*’*U 

e.  Go  back  to  step  II. 

9.  For  S  =  1  to  5ta(al : 

Update  the  Lagrange  parameters:  =  XI  * 

10.  Go  back  to  to  step  II.B 


Algorithm  2:  Description  of  the  dynamic  augmented  Lagrangian  method. 
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The  parameters  in  Algorithm  2  are  defined  by: 


in  the  time  domain  by  the  Neimiark  method.  Note  that  for  linear  elastic  prob¬ 
lems,  steps  which  involve  forming  could  be  moved  out  of  the 

loops  in  both  Algorithms  1  and  2. 


Below,  the  two  spring  problem  is  again  considered  for  a  suddenly  applied 
load.  A  comparison  with  the  exact  results  in  Figure  10  shows  excellent  agree¬ 
ment. 


2 

1 

2k  +  fi  -(k  +  fi) 

3^ 

6^ 

-(k  +fi)  fc  +  /X 

1 

6^ 

1 

3^1 

This  system  leads  to  the  following  collection  of  matrices: 
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The  initial  conditions  are  assumed  to  be: 


0^ 

.8 

Oti, 

7 

“"2 

“  [oj 

Oli, 

*  [0] 

OiiB 

-2.4 

7 

Furthermore,  the  system  constants  are  taken  to  be: 


(112) 


(113) 


(114) 


(115) 


a 


*  =  1 


TJl  =  1 


q  =  -.1 


(116) 
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Figure  10:  Comparison  of  exact  and  augmented  Lagrangian  results. 


4.  Examples 

This  section  contains  examples  of  the  augmented  Lagrangian  solution 
method  for  both  the  static  and  dynamic  cases.  The  examples  presented  here 
are  compared  to  work  done  by  previous  researchers.  In  particular,  the  aug¬ 
mented  Lagrangian  method  is  compared  to  the  results  of  references  [6]  and  [7] 
for  the  perturbed  Lagrangian  approach,  reference  [15]  for  the  dynamic  case 
and  reference  [14]  for  the  Hertz  solutions. 

The  augmented  Lagrangian  method  was  applied  to  the  system  of  a  rigid 
punch  into  an  elastic  foundation,  illustrated  in  Figure  11.  The  finite  element 
mesh  is  successively  refined  in  order  to  get  an  idea  on  the  convergence  charac¬ 
teristics  of  the  system.  Using  information  from  reference  [7],  the  following 
parameters  were  used; 

^funeh  ~  •'pufieA  “  0.0  E/mnpaften  —  10®  /ownAo^Wn  ”  0.3  (117. a) 

l*»u»eA  “  20  It  ^  foundatim  ~  22  (ll7.b) 


Finite  element  model  Displaced  shape 


Log  (  Penalty  Value  ) 


Figure  13:  Comparison  of  converged  results  for  the  3  element  mesh. 

Figures  12  and  13  display  the  results  of  applying  the  penalty,  perturbed 
Lagrangian  and  augmented  Lagrangian  methods  to  the  three  block  problem. 
Figure  12  shows  the  results  after  1  iteration  of  each  of  the  methods.  Note  that 
the  augmented  Lagrangian  method  gives  accurate  results  over  a  much  wider 
range  of  penalty  parameter  values.  The  converged  results  are  shown  in  Figure 
13.  In  this  case,  the  augmented  Lagrangian  approach  always  }rields  the  correct 
results,  while  the  other  methods  only  approximate  the  solution  for  smaller 
penalty  values.  Table  2  displays  the  comparison  on  the  number  of  iterations  for 
the  three  methods,  using  a  penalty  value  of  fi-  10'^.  Observe  that  even  though 
the  solutions  are  in^>roved  over  the  standard  penalty  and  perturbed  Lagran¬ 
gian  methods,  the  values  in  Table  2  show  that  little  extra  work  is  required. 


Method 

Mesh  1 

penalty 

3 

perturbed  Lagrangian 

2 

augmented  Lagrangian 

2 

Mesh  2  Mesh  3 
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Figure  14:  Finite  element  mesh  for  the  Hertz  contact  problem. 

In  order  to  further  illustrate  the  augmented  Lagrangian  method,  the  finite 
element  system  in  Figure  14  was  studied.  This  mesh  represents  a  sphere  in  con¬ 
tact  with  a  rigid  wall.  The  parameters  for  this  system  are: 

E=108  „s:0.3  /»  =  0.01  R  =  B.0  (118) 

The  study  of  a  static  analysis  involved  three  different  load  values  applied  to  the 
top  of  the  sphere.  Results  of  the  loading,  compared  to  the  classical  Hertz  solu¬ 
tion  are  illustrated  in  Figure  IS.  An  iteration  count  comparing  the  three 
methods  can  be  found  in  Table  3.  Each  of  these  method  used  a  penalty  value  of 
/i  =  10'^.  A  dynamic  analysis  was  also  carried  out.  assuming  that  the  sphere  was 
given  an  initial  velocity  of  0.1  downward.  The  results  displayed  in  Figure  16 
represent  the  maximum  impact  values.  In  this  dynamic  analysis,  only  the  aug¬ 
mented  Lagrangian  method  obtained  reasonable  results. 
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5 
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Figure  15;  Static  results  for  the  Hertz  contact  problem 


S.  Recommendations  for  further  work 

This  section  contains  a  list  of  possible  extensions  to  the  solution  method 
and  its  associated  algorithms.  In  particular,  proposals  are  made  which  would 
improve  the  rate  of  convergence  of  the  algorithm,  make  the  slide  line  computa¬ 
tions  more  efficient  and  extend  the  algorithm  to  cover  new  problem  domains. 

Using  equation  (93),  a  second  order  updating  scheme  could  be  added  to 
the  augmented  Lagrangian  method.  By  comparing  equations  (83)  and  (84),  the 
following  scheme  could  be  developed. 

-(9«d(A*)]"‘  =  V®p(u)  +  fjLl  (119) 

Thus,  the  iteration  corresponding  to  equation  (93)  is  obviously; 

X*  ^  ,  =  Xfc  +  (u)  +  Mljvd  (X* )  (120) 

As  the  computation  proceeds,  values  are  found  for  p  (u)  and  Vp  (u)  =  X*  +  /iu. 
These  values  could  then  be  used  to  generate  an  approximation  to  V^p(u).  The 
need  to  use  second  derivatives  is  eliminated  entirely  if  one  chooses  to  solve  the 
contact  problem  by  some  quasi-Newton  methods  such  as  the  DFP  or  BFGS.  See 
reference  [4]  for  more  details  on  these  methods. 

A  further  enhancement  might  be  to  extend  the  slide  line  contact  logic.  As 
was  mentioned  in  Section  3,  the  logic  is  not  well  defined  when  multiple  segments 
are  in  contact  with  a  single  surface.  For  four  node  isoparametric  elements,  the 
extension  is  straightforward.  However,  an  extension  to  cover  three  dimensional 
problems  would  entail  the  use  of  a  clipping  algorithm  to  determine  if  segments 
were  in  contact.  Essentially,  the  algorithm  is  cubic  in  the  number  of  nodes  for 
this  case,  but  if  linear  elements  are  used,  a  proper  implementation  could 
reduce  the  required  computation  significantly.  The  extension  to  cover  higher 
order  elements  is  less  clear.  It  is  believed,  however,  that  by  using  techniques 
from  constructive  solid  geometry  and  exploiting  the  convex  hull  properties  of 
the  finite  element  mesh,  a  reasonably  efficient  method  could  be  developed. 

A  third  area  which  needs  further  work'  is  the  development  of  algorithms 
which  include  friction  on  the  coritact  surface.  In  particular,  methods  which 
develop  consistent  tangent  matrices  need  to  be  devised.  The  addition  of  fric¬ 
tional  contact  to  the  algorithms  presented  here  would  greatly  enhance  their 
generality. 

Finally,  the  contact  algorithms  should  not  be  limited  to  small  deformation 
kinematics.  All  of  the  solution  efforts  that  have  been  described  here  need  to  be 
evaluated  for  situations  in  which  motions  on  the  contact  surface  are  finite. 
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6.  Summary  and  Conclusions 

An  augmented  Lagrangian  method  for  the  solution  finite  element  contact 
problems  offers  several  advantages  over  the  standard  penalty,  classical  Lagran¬ 
gian  or  perturbed  Lagrangian  methods.  By  using  the  augmented  Lagrangian 
method,  the  user  has  direct  control  over  the  accuracy  of  a  given  problem,  since 
it  is  the  number  of  iterations  that  control  the  tolerance  of  the  solutions.  The 
value  of  the  penalty  parameter  controls  the  rate  of  convergence.  Expressions 
have  been  presented,  based  on  error  considerations,  which  may  be  used  as  a 
guide  to  selecting  a  reasonable  penalty  parameter. 

The  algorithm  for  the  augmented  Lagrangian  method  can  be  simplified  to 
incorporate  both  the  standard  penalty  and  perturbed  Lagrangian  Methods. 
Furthermore,  this  method  is  applicable  to  both  static  and  dynamic  problems. 
The  results  presented  in  section  4  show  that  the  method  compares  very  well  to 
work  done  by  previous  researchers.  Convergence  is  quite  fast,  even  in  cases 
where  the  other  methods  fail. 
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